function th = applyfdr(stats, stype, levels, df1, df2, lnmeth)
% applyfdr  - apply FDR thresholding to given statistic
%
% FORMAT:       th = applyfdr(stats, stype, levels, df1 [, df2 [, lnmeth]])
%
% Input fields:
%
%       stats       statistical values (p/r/t/F)
%       stype       either of {'p'}, 'r', 't', 'F'
%       levels      q(FDR) levels
%       df1         d.f. 1
%       df2         d.f. 2
%       lnmeth      if given and true (boolean!), use
%                   c(V) = ln(V) + E  instead of  c(V) = 1 approximation
%
% Output fields:
%
%       th          thresholds in original statistic
%
% See also fdr_thresholds

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
   (~isa(stats, 'double') && ...
    ~isa(stats, 'single')) || ...
   ~ischar(stype) || ...
    numel(stype) ~= 1 || ...
   ~any('fprt' == lower(stype)) || ...
   (~isa(levels, 'double') && ...
    ~isa(levels, 'single')) || ...
    isempty(levels) || ...
    any(isinf(levels(:)) | isnan(levels(:)) | levels(:) <= 0 | levels(:) >= 1) || ...
   (lower(stype) ~= 'p' && ...
    ((~isa(df1, 'double') && ...
      ~isa(df1, 'single')) || ...
     numel(df1) ~= 1 || ...
     isnan(df1) || ...
     df1 < 1 || ...
     (lower(stype) ~= 'f' && ...
      df1 < 2) || ...
     (lower(stype) == 'f' && ...
      (nargin < 5 || ...
       (~isa(df2, 'double') && ...
        ~isa(df2, 'single')) || ...
       numel(df2) ~= 1 || ...
       isnan(df2) || ...
       df2 < 1))))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 6 || ...
   (~islogical(lnmeth) && ...
    ~isnumeric(lnmeth)) || ...
    numel(lnmeth) ~= 1
    lnmeth = 0;
else
    lnmeth = double(lnmeth);
end
nstats = numel(stats);
levels = levels(:);
stats = stats(:);
stats(isnan(stats) | stats == 0) = [];

% for empty stats, set to Bonferroni level - eps
if isempty(stats)
    
    blevels = (1 / nstats) * levels - eps;
    
    % depends on type
    switch (lower(stype))
        case {'f'}
            th = custom_finv(1 - blevels, df1, df2);
        case {'p'}
            th = blevels;
        case {'r'}
            th = correlinvtstat( ...
                custom_tinv(1 - 0.5 * blevels, df1), df1 + 2);
        case {'t'}
            th = custom_tinv(1 - 0.5 * blevels, df1);
    end
    if lnmeth
        th = th(:, [1, 1]);
    end
    return;
end

% depends on type
switch (lower(stype))
        
    % F stats
    case {'f'}
        th = custom_finv(1 - fdr_thresholds( ...
            1 - custom_fcdf(abs(double(stats(:))), df1, df2), ...
            levels(:), lnmeth), df1, df2);
    
    % already p values
    case {'p'}
        th = fdr_thresholds(stats(:), levels(:), lnmeth);
        
    % correlation r
    case {'r'}
        th = correlinvtstat(custom_tinv(1 - fdr_thresholds( ...
            1 - custom_tcdf(abs(correltstat(stats(:), df1 + 2)), df1), ...
            levels(:) / 2, lnmeth), df1), df1 + 2);
        
    % t stats
    case {'t'}
        th = custom_tinv(1 - fdr_thresholds( ...
            1 - custom_tcdf(abs(double(stats(:))), df1), ...
            levels(:) / 2, lnmeth), df1);
end
